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ABSTRACT 

Cosmological constraints from cluster surveys rely on accurate mass estimates from the mass- 
observable relations. In order to avoid systematic biases and reduce uncertainties, we study the 
form and physical origin of the intrinsic scatter about the mean Sunyaev-Zel'dovich (SZ) flux-mass 
relation using a hydrodynamical simulation of galaxy cluster formation. We examine the assumption 
of lognormal scatter and detect non- negligible positive skewness and kurtosis (> 0.5) for a wide range 
of limiting masses and redshifts. These higher-order moments should be included in the paramctriza- 
tion of scatter in order not to bias cosmological constraints. We investigate the sources of the scatter 
by correlating it with measures of cluster morphology, halo concentration, and dynamical state, and 
we quantify the individual contribution from each source. We find that statistically the impact of 
dynamical state is weak, so the selection bias due to mergers is negligible. On the other hand, there 
is a strong correlation between the scatter and halo concentration, which can be used to reduce the 
scatter significantly (from 12.07% to 7.34% or by ~ 40% for clusters at z = 0). We also show that a 
cross-calibration by combining information from X-ray followups can be used to reduce the scatter in 
the flux-mass relation and also identify outliers in both X-ray and SZ cluster surveys. 
Subject headings: dark matter — galaxies: clusters: general — hydrodynamics — intergalactic medium 
— methods: numerical 



1. INTRODUCTION 

The evolution of structure in the Universe is thought to 
be a hierarchical process driven by gravitational instabil- 
ity acting on primordial density fluctuations. In this dy- 
namical process, smaller clumps of matter merge to form 
bigger ones within a "cosmic web" of spatial structure in- 
corporating matter within tenuous sheets, higher-density 
filaments, and large matter concentrations at the nodal 
points of the web. Clusters of galaxies occupy the top of 
the mass hierarchy, being the largest objects that have 
had time to collapse and form due to their self-gravity. 
Therefore, they are the most recent structures to form. 
Clusters have two very important features: they are very 
massive (1O 14 ~ 15 M0) and populate the exponential tail 
of the mass function , which is sensitive to dark energy 
(|Haiman et al.|[200ll ). and because of their deep potential 
wells, they are "matter traps," preventing their internal 
constituents from escaping. The first aspect, combined 
with their young age, makes clusters an excellent probe of 
cosmology, being especially sensitive to dark matter and 
dark energy. The second makes clusters excellent labora- 
tories for studying key astrophysics problems such as star 
and galaxy formation. Because of their central impor- 
tance to both cosmology and astrophysics, clusters are 
being targeted by a variety of observational programs. 

Clusters can be detected optically via starlight from 
their member galaxies, in X-ray due to thermal emission 
from hot ga s, and using the thermal Su nyaev-Zcrdovich 
(SZ) effect (ISunvaev fc Zeldovichl [19721) . The SZ effect 
arises from the inverse- Compton upscattering of cos- 
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mic microwave background (CMB) photons by the hot 
electrons in a cluster (several keV), leading to an in- 
crease of the CMB brightness temperature at frequencies 
above about 250 GHz and a decrease at lower frequen- 
cies. This distortion of the Planck spectrum can be mea- 
sured by se nsitive, high-resolution CMB telescopes from 
the ground dHincks et alj|2009t [Staniszewski et alj|2009t 
iVanderlinde et al.l I2010T ) . and it is one of the best way s 
to find clusters at higher redshifts (jBrodwin et al.l[20~lCf ) , 
since the SZ effect is (almost) redshift-independent. 

Upcoming surveys will detect clusters using their SZ 
signatures. However, instead of mass M, SZ surveys 
measure Y, the SZ Compton optical depth (the so-called 
y-decrement) integrated over a portion of a cluster's pro- 
jected area. To obtain the cluster mass distribution 
from surveys, we need a so-called mass-observable (Y- 
M) relation. Using the mass distribution to constrain 
cosmological parameters requires that we know the er- 
rors introduced in this process, so we must also quan- 
tify the scatter and the bias in the mass-observable re- 
lation as functions of redshift (|Lima !Hlul[200l [2005; 
ICunha fc Evra7d1[20l0l) . 

Because of the exponential shape of the cluster mass 
function, scatter in the X-M relation (for any observable 
X that is positively correlated with M) boosts the num- 
ber density of clusters observed in logarithmic bins of X, 
as the overall number of lower-mass clusters scattering to 
higher values of X far exceeds the number of high-mass 
cl usters scattering in the opp osite direction. As shown 
in iBhattacharva et al.l ()2010[ ). a few percent systematic 
error in mass, which can arise due to bias in the X-M 
relation, leads to a significant difference in the tail of the 
mass function. Thus, misestimating the scatter and bias 
in the X-M relation c an lead to biases in the cosmo- 
logica l constraints (e.g. iRandall et al.l 120021 : iWik et al~l 
120081 ). Moreover, previous work forecasting cosmologi- 
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cal constraints based on cluster survey s has assumed the 
distribution of scatter to be lognormal (|Lima "fcHul[2005l: 
Cunha 200jJ. If this assumptio n were invalid, it w ould 
also lead to biases in the results (jShaw et alll2010"al ). 

Understanding the physical sources of scatter can help 
us to reduce this scatter and improve the mass estimates. 
For example, the use of core-excised quantities reduces 
the effects of cool cores in clusters and hence also the 
large scatter i n the X-ray luminosit y-temperature [Lx~ 
T x ) relation (I Allen & F^bianl 119981 ) or the L x -M rela- 
tion (jMantz et all 120101 ) . Using the halo concentration 
as a third p arameter can red uce the scatter in the Tx~ 
M relation (|Yang et al.l 120091) . Combinations of observ- 
ables with oppositely trending scatter can also help, as 
shown by the tight correlation between the X-ray coun- 
tcrpar t of the Comptqn ^-pa rameter, Yx = M gas Tx, and 
mass (jKravtsov et al.1 [200(3 ). These examples illustrate 
the possibility of obtaining better mass estimates if our 
knowledge of the physical origin of scatter can be im- 
proved. 

Because of spatial resolution constraints, the Y—M 
relation has been studied using two types of hydro- 
dynamic simulations. One type includes only adia- 
batic physics but has sufficient statistics to quantify 
the scatter and the b ias in the mass-observable relation 
(|Hallman et al.l 120071 ). The other includes extra physics 
such as ra diative cooling, sta r formation, and supernova 
feedback (iNaeai et al.l 120071) or feedback from quasars 
(jBhattacharva et all 120081 ). but with limited statistics 
(16 and 10 halos, respectively, in the cited references). 
Another approach is to include gas physics using a semi- 
analytic gas prescription in halos obta ined from a dark- 
matt e r only (DMO) sim ulation (e.g. iBode et all 120071 
l2009t IShaw et al.l l2010bh . In particular, |Shaw_et3l 
(|2008l I2010a| ) have quantified the Y-M relation and 
its non-gaussianity using th e "DMO+semiana lytic" ap- 
proach. As pointed out by IShaw et all ()2008f ) , the hy- 
drodynamic simulations tend to show slightly different 
scatter compared to the "DMO+semianalytic" case es- 
pecially for overdensitics ^500 times the critical density 
of the universe. These differences need to be understood. 
Only recently have there been attempts to incorporate 
extra baryonic physics into a large cosmological simula- 
tion to study the in trinsic variances in t he scaling rela- 
tions. For instance, iStanek et al.l (|2010l ) have used res- 
imulations from the Millennium Gas Simulation with gas 
dynamics treated in both gravity-only and cooling plus 
preheating prescriptions to study the scaling relations 
and correlations among cluster structural properties and 
observables in X-ray and SZ. 

In this study we begin a systematic investigation of 
the physical sources of the intrinsic Y-M scatter and 
their impact on the form of scatter for the purpose of 
improving our knowledge of both cluster formation and 
cluster cosmology. We use cosmological simulations to 
obtain clusters with sufficient statistics and to incorpo- 
rate hydrodynamic processes including interactions with 
large-scale environment, variations in cluster structures, 
and merger-induced shock heating and departure from 
hydrostatic equilibrium. Since we would like to focus 
on the scatter driven by gravitational effects only, radia- 
tive cooling and heating mechanisms are not included. 
We will address the influence of radiative cooling and 
feedback explicitly in a separate paper. In this work, 



we examine the assumption of lognormal scatter, which 
has important implications for self-calibration studies of 
cluster surveys. We investigate various sources of scatter, 
such as halo concentrations, dynamical state, and cluster 
morphology, by correlating the scatter with quantitative 
measures of each source. We show that these correlations 
can be used to reduce the scatter and tighten the scaling 
relation. We also discuss possible applications and issues 
when combining SZ and X-ray scaling relations. 

The outline of this paper is as follows. In § [2] we sum- 
marize the key components of our simulation and anal- 
yses, including the numerical methods and simulation 
parameters, a brief overview of the SZ effect, merger tree 
construction, and how we create idealized cluster samples 
to disentangle the sources of scatter. The Y-M relation 
and the form of its scatter arc presented in § [3] Possible 
sources of scatter are investigated in § [4j In § [5] we ex- 
plore the possibility of combining SZ and X-ray scaling 
relations to improve cluster mass estimates. Finally, we 
discuss our results and give the conclusions in § [6] 

2. METHOD 

2.1. Simulation 

The simulation described here was performed using 
FLASH, an Eulerian hydrodynamics plus iV-body code 
which has been applied to a wide range of problems and 
exten sively validated for hvdrod vnamical (|Calder et all 
2002 ) and cosmological A-body (|Heitmann et alll2005l 
20081 ) applications. We used version 2.4 of FLASH to- 
gether with the loc al trans f orm-b ased multigrid Poisson 
solver described by iRickerl (120081 ) . Because we are con- 
cerned in this paper only with the effect of gravity-driven 
processes on mass-observable relations, the calculation 
described here did not employ radiative cooling or feed- 
back due to star formation or active galaxi es. Here we 
jve a brief summary and refer the readers to lYang et all 
20091) for details of the numerical methods and merger 
tree analysis. 

The results presented here are based on a FLASH 
simulation of structure formation in the ACDM cosmol- 
ogy within a 3D cubical volume spanning 256ft. -1 Mpc. 
Initial conditions were gener ated for a starting red- 
shift z = 66 using GRAFIC (jBertschingerl 120011 ) with 
an initial power spectrum generated using CMBFAST 
(jSeliak fc Zaldarriagalll996l ). The cosmological param- 
eter values used were cho sen to be consistent with the 
third-year WMAP results (jSpergel et alll2007f l: present- 
day matter density parameter fl m Q — 0.262, present-day 
baryonic density parameter Qto = 0.0437, present-day 
cosmological constant density parameter J7ao = 0.738, 
matter power spectrum normalization erg = 0.74, and 
Hubble constant h = 0.708 (H = 100ft km s" 1 Mpc -1 ). 
The simulation contains 1024 3 dark matter particles with 
a particle mass m p — 9.2 x lO 8 /i _1 Af . The mesh used 
for the gasdynamics and potential solution was fully re- 
fined to 1024 3 zones, which corresponds to a comoving 
zone spacing of 250/i _1 kpc. Considering the effect of 
resolution on the computed ab undances of halos of dif- 
ferent mass (|Lukic et al.l 120071) . with these parameters 
we are able to capture all halos containing more than 
3150 particles (i.e. total mass 2.9 x 1O 12 /i _1 M ) and 1150 
particles (i.e. 1.1 x l0 12 h~ 1 M G ) at z = and z = 1, 
respectively. The halos arc identified using the friends- 
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of-friends (FOF) algorithm. The overdensity mass and 
radius, Ma and Ra, are then found by growing spheres 
around each FOF center until the averaged total density 
is A times the critical density of the universe. The sim- 
ulation was carried out using 800 processors of the Cray 
XT4 system at Oak Ridge National Laboratory, requir- 
ing a total of 16,500 CPU-hours. 

2.2. Sunyaev-ZcVdovich Effect 

The thermal Sunyaev-Zel'dovich effect is caused by in- 
verse Compton scattering of CMB photons off the hot 
electrons inside galaxy clusters. Assuming that rela- 
tivistic corrections are small, the resulting distortion of 
the CMB temperature can be written as AT/Tcmb = 
g v (x)y where g v (x) = x(coth(x/2) — 4) with x = 
hv/ksTcMB, ks is the Boltzmann constant, v is the 
frequency of observation, and Tqmb is the mean CMB 
temperature at the current epoch. Here y is the Comp- 
ton y-parameter, which can be written as 



y 



m e c 2 



n e {l)T e (l)dl, 



(1) 



where n e (l) is the electron density profile, ot is the 
Thomson scattering cross section, and the integration 
is along the line of sight, which is defined to be the x- 
direction in the simulation box. 

Note, however, that the observable is the integral of 
the temperature distortion over the cluster's projection 
onto the sky. For a cluster at redshift z, it is given by 



Y(M,z) 



1 k,B<JT 

dA(z) 2 m e c 2 



n e (l)T e (l)dV, (2) 



where dA{z) is the angular diameter distance to the clus- 
ter, and the integration is over the volume of the cluster. 
In the literature, sometimes the factor is omitted 

from Eq. [2J In this study, when we omit the factor l/d A , 
we denote the temperature distortion as Y(M); other- 
wise it is denoted as Y(M,z). Note that Y(M,z) is 
dimensionless, while the units of Y(M) are Mpc 2 . In 
the following sections in which we investigate the distri- 
bution and origin of intrinsic scatter, we adopt Y(M) 
unless explicitly stated otherwise. 

Following the convention in the literature, we denote 
Ya as the SZ flux integrated out to certain overdensity 
radius Ra- In the simulation box, cell-averaged informa- 
tion about the gas density and temperature is stored for 
each grid cell. Thus for each cluster the integrated SZ 
distortion is calculated using 



Y A (M,z) 



1 kB<JT 



d A {z) 2 m e c 



r p <R& 



n e! iT eti AVi, (3) 



where the summation is over all the grid cells across the 
cluster volume within a cylinder of projected radius Ra, 
and AVi is the volume of each grid cell. 

2.3. Merger Tree Analysis 

In order to explore the influence of cluster formation 
and merger history on the intrinsic scatter, we construct 
a merger tree for each cluster in our simulation in the 
following way. Our simulation generates snapshots that 
contain particle tags and positions every 100/i -1 Myr be- 
ginning at z = 2. For each snapshot, all the groups with 



more than 10 particles are found using the FOF halo 
finder with linking length parameter b = 0.2. Between 
successive outputs at times t = t n , we find the progen- 
itors at time i„_i for all the halos at t n by tracing the 
particle tags, which are uniquely assigned to each particle 
at the beginning of the simulation. For each halo at t n 
we record the masses of its progenitors, the masses they 
contributed to the halo, and the number of unbound par- 
ticles. Then the merger trees are constructed by linking 
all the progenitors identified in the previous outputs for 
halos above our halo completeness limit at z = 0. De- 
riving the mass accretion histories is straightforwardly 
accomplished by following the mass of the most mas- 
sive progenitor back in time. Cluster formation time is 
often defined as the epoch when a cluster exceeds a cer- 
tain fraction of its final mass. The commonly-adopted 
thresholds include 10%, 25%, 50%, and 70%. In discus- 
sion below we present results using the 50% threshold. 

To directly quantify the dynamical state of clusters 
without relying on morphology, we find the time since 
last merger for each cluster in our simulation. In our 
analysis, mergers arc defined in two ways: the mass-jump 
definition, in which a merger is present if there is a mass 
jump larger than some threshold in the halo's assembly 
history; and the mass-ratio definition, which identifies a 
merger if the ratio of contributed masses from the first- 
and second-ranked pr ogenitors is less than a certain value 
(jCohn fc W hite 2005j) . To study the variations in cluster 
observables induced by different types of mergers, we use 
1.2 and 1.33 as thresholds for the mass-jump definition 
and 10:1, 5:1, and 3:1 in the mass-ratio definition. In 
this paper, by 'merging clusters' at a given lookback time 
we will refer to those identified by at least one of these 
five merger diagnostics in the preceding 3 Gyr, chosen to 
be long enough such that mergers with different impact 
parameters and mass rat ios would have ret urned to virial 
equilibrium within -R500 (jPoole et al.ll2006T ) . The mergers 
are 'major' if the mass jump is larger than 1.2 or if the 
mass ratio is less than 5:1; 'minor' mergers, on the other 
hand, have mass ratios between 10:1 and 5:1. 

2.4. Idealized Cluster Samples 

One of our main goals is to investigate the possi- 
ble sources of intrinsic scatter, including the variations 
due to concentration, departure from hydrostatic equi- 
librium, merger boosts, and cluster morphology. In or- 
der to distinguish the contribution from each source, we 
construct a set of idealized cluster samples with different 
assumptions. Starting from a sample with the most pos- 
sible constraints, we add one source of scatter at a time. 
By comparing the scatter of the idealized samples and 
the simulated sample, we can tell whether the scatter 
can be successfully reconstructed, or yet other sources 
still need to be found. 

To this end we construct five idealized samples, going 
from sample A, with the most constraints, to sample E, 
which includes the most sources of variation. For sample 
A, only the cluster masses are taken from the simulation. 
Given the mass, the halo concentra tion c is compute d 
using the best- fit c-M relation from IShaw et al.l ((2006ft . 
With the mass and halo concentration, the total den- 
sity and gas density are as signed using the NFW profile 
(Navarro, Frenk, & White 119951: 119961 hereafter NF W) 
and a corc-softcncd NFW profile (|Siiacki et al.ll2007l ) , re- 
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TABLE 1 

Summary of models used for constructing the idealized cluster samples. 



Sample 


Assumptions / Constraints 


Sources of Variation 51 


A 


Spherical + f g + HSE^_ + c(M) 


None 


B 


Spherical + f g +"HSE 


c 


G 


Spherical + f g + No Merger Boost 


c + Random Gas Motion 


D 


Spherical + f g 


c + Random Gas Motion + Merger Boost 


E 


Spherical 


c + Random Gas Motion + Merger Boost + f g 


S 


Simulated 


All 


SS 


Simulated + Spherical 


All - Morphology 



All include scatter resulted from particle shot noise, finite cell resolution, and the simulated observation procedure. 
Q Hydrostatic equilibrium enforces no merger boost and no random gas motion. 



spectively. In the core-softened NFW profile, the core 
radius is set to be 0.02 times the virial radius, and the 
gas fraction, f g = 0.12, is also fixed. The pressure and 
temperature profiles arc then computed assuming hydro- 
static equilibrium (HSE). Finally, dark matter particles 
and zone-averaged cell quantities are assigned according 
to the profiles, assuming spherical symmetry. They are 
stored using the same file format as the simulated clus- 
ters, allowing them to be analyzed in the same way as 
the simulated clusters. Any scatter in the resulting Y—M 
relation can only be due to particle shot noise, finite cell 
resolution, and the simulated observation procedure. 

Sample B is generated using a similar procedure, ex- 
cept that the c-M relation is assumed to have a log- 
normal distribution with a dispersion of 0.22 (| Jindl200Clt 
IBullock et al"][200ll IDolae et alJl200i IShaw et al J 1 2006ft . 
In other words, the variation in concentration should be 
the only additional source of Y-M scatter for sample B. 

In addition to the variation in concentration, clusters 
in sample C are allowed to depart from HSE due to pres- 
sure support from random gas motions after mergers. For 
each cluster, we compute the gas velocity dispersion's 
radial profile from the simulation and include an extra 
pressure term, P ran d = per 2 , in the HSE equation for 
computing the thermal pressure. In this way, we are ef- 
fectively taking into account the incomplete virialization 
after merger events as another source of scatter. 

With sample D, we model the effect of mergers by 
including not only incomplete relaxa tion but also the 
merg er boosts due to shock heating (jRicker fc Sarazinl 
120011 iPoole et al .112007ft . Therefore, we extract the time 
histories of mass and integrated SZ flux during mergers 
from the ideal merger simulations in iPoole et al.l (|2007ft 
and apply a boost in Y to each idealized cluster using 
the actual times since last merger in the simulation and 
the mass ratios of th ose mergers. Note t hat the time evo- 
lution of mass from IPoole et al.l (|2007ft is measured for 
overdensity A = 500, but the integrated SZ flux is only 
available for A = 2500. Thus the Y-M scatter of sample 
D should be considered as an upper limit to the effects 
of merger boosts. 

Sample E adds variation in the gas fraction by using 
the actual gas fraction computed for each simulated clus- 
ter. Sample E essentially includes all possible sources 
of scatter investigated in this paper except the effect of 
cluster morphology. In § !4.5l we will compare results from 
these idealized samples to the simulated clusters (Sam- 
ple S). Because the idealized samples are all constructed 
under the assumption of spherical symmetry, we derived 
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Fig. 1. — Normalization and slope of the Y— M relation as func- 
tions of redshift. The dashed lines are the self-similar prediction. 

another sample (Sample SS) using gas profiles directly 
extracted from the simulated clusters, such that it in- 
cludes all sources of scatter except the morphological ef- 
fect. These models are summarized in Table [1] 

3. THE Y-M SCALING RELATION 

3.1. Normalization and Slope 

For each simulation output between z = and z = 
1.5, we derive the Ya — Ma relation for all clusters with 
M 500 > 2 x 10 13 Af o (619 clusters at z = and 223 
clusters at z = 1) and fit it with a power law of the form 

where the normalization at 10 14 /i _1 AfQ in units of 10~ 6 , 
An, and slope a are found by fitting the data points 
using a Levenberg-Marquardt algorithm. An overdensity 
of A = 500 is adopted throughout the paper. 

The normalization (relative to the value at z = 0), 
slope, and scatter as functions of redshift arc plotted in 
Figure [U As also found in previous adi abatic simulations 
(jda Silva et al.l 120041: iMotl et al.l [2005ft . the evolution of 
the normalization and slope is consistent with the self- 
similar prediction, a = 5/3 and Ai^z) oc E(z) 2 / 3 , where 
E(z) = [O m o(l + Z ) A + ^ao] 1 ^ 2 , though we cannot rule 
out the case of no evolution because the small number of 
high-mass clusters at higher redshift limits the constrain- 
ing power of the data. We find ^14(2 = 0) = 5.43 ±0.47, 
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Fig. 2. — RMS scatter as a function of Msoo- Different curves 
arc for different rcdshifts. In general, the scatter decreases with 
both mass and redshift. See Eq.[7]for the best-fit relation. 



which is in agreement with the adiabatic ru ns in | Nagai 
2001) (A u (z = 0) = 4.99, / fc = 0.14) and IShaw et al 



20081 ) (A 14 (z = 0) = 4.07, f b = 0.11) after taking into 



account the differences in the baryon fraction used in the 
simulations (fb = 0.167 in our simulation). 

3.2. Scatter in the Y-M Relation 

The RMS scatter around the best-fit relation is defined 
for N clusters as 



o~ym 



^(log^-log^) 2 



N — 1 



1/2 



(5) 



where Yi is the measured flux of the ith cluster, and Yi is 
the flux predicted by the best-fit relation for that cluster. 
Hereafter, we use the notation 



5 log Y = log Y — log Y 



(6) 



for the deviation from the mean relation for each cluster. 
For each redshift from z = to z = 1, we compute 
the RMS scatter for 5 mass bins and plot the result in 
Figure [2] (the data for some of the redshifts are omitted 
for clarity). The scatt er is ^ 5 — 15%, consistent with 
previous findings (e.g. [Nagai 2006T). Moreover, we find 
that in general the scatter decreases with both mass and 
redshift. We fit the scatter using the functional form 



a(M,z) = AlogM + Blog(l + z) + C, 



(7) 



where the best-fit coefficients are A = —7.06 ± 0.28, 
B = -11.20 ± 0.81, and C = 7.70 ± 0.19. The mass de- 
pendence may be due to increasing non-lognormality of 
the scatter when considering low-mass clusters (see § 13.31 
for details). The redshift evolution may be understood 
by considering the self-similar model, in which all quanti- 
ties for collapsed objects can be expressed in terms of the 
characteristic mass scale, M* oc (1 + z)~ 6 /(™+ 3 ), where 
n is the spectral index of the scale-fr ee primordial power 
spectrum, P(k) oc fc™ (jKaiserl [198(1 . Therefore, Eq. [7] 
is equivalent to the expression a = A' log(M/M*) + B' . 
Note that Eq. [7] is the first attempt in the literature to 
quantify the scatter using a functional form of mass and 
redshift. This expression should be useful for future stud- 
ies that require assumptions about the form of scatter. 



3.3. Non-Lognormal Scatter 

Since the SZ flux is only proportional to the first power 
of gas density, it is sensitive to the contribution from low- 
density gas clumped along the line of sight to but not 
associated with a cluster. This is in contrast to the X- 
ray luminosity, which is proportional to the square of gas 
density. For this reason the Y-M relation has been found 
to h ave a high-scatter t ail in the distributio n of its scat- 
ter ijWhite et all 120021: iHallman et ai1l2007| ). Since any 
deviations from the lognormal scatter would bias cosmo- 
logical constraints based on cluster counts (|Shaw et al.l 
2010a). we would like to examine whether the form of the 
log scatter can be well approximated by a Gaussian dis- 
tribution, or whether generalizations of the parametriza- 
tion need to be considered. 

A purely Gaussian distribution can be described ex- 
actly using only its mean value [i and variance a : 



G(x) 



1 



cxp 



(x - /i) 5 



2a 2 



(8) 



A nearly-Gaussian distribution 
imated using the Edgew or th 
Bernardeau fc Kofmanl 119951 ; iBlinnikov fc Moessnerl 
19981) . 



be approx- 
expansion (e.g. 



G(x) « G(x) 



7 d 3 G k d*G 7 2 d 6 G 
~ 6 ~dx T + 24 dx 7 + 72 "^6"' 



(9) 



which is parametrized by four moments - the mean and 
the variance describing the Gaussian distribution, plus 
the skewness (7) and the kurtosis (k) describing the de- 
viation from gaussianity. The skewness is defined as 



7 



((z-M) 3 } 



and the kurtosis as 



K = ; 3. 



(10) 



(11) 



We compute the skewness and kurtosis of the Y-M 
scatter for clusters at z = above different mass thresh- 
olds and plot the results in Figure |31 The error bars 
represent the uncertainty due to finite sample size and 
are estimated using 10 3 Monte-Carlo realizations of ran- 
dom sampling from a nearly-Gaussian distribution given 
the measured skewness and kurtosis as in Eq. [S] Note 
that because of finite sample size the measured skewness 
and kurtosis would underestimate the intrinsic values of 
the underlying distribution. This bias is represented by 
the offset between the data points and the middle points 
of the error bars. 

We find that the scatter is non-lognormal with posi- 
tive skewness and kurtosis when including only massive 
clusters (Msoojim > 10 14 M Q ) or when more and more 
low-mass clusters are included (M 50 o,ii m <5x 10 13 M Q ). 
Due to the limited number of clusters at the high-mass 
end, the log scatter there is expected to follow Poisson 
statistics and deviate from a Gaussian form. For the 
lower mass range, on the other hand, we find that there 
is a tail toward positive values in the distribution of scat- 
ter, which increases the skewness and kurtosis. By visual 
inspection of clusters in the tail, we find that these ob- 
jects happen to have elongated shapes or clumped gas 
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Fig. 3. — Skcwncss and kurtosis of the Y-M scatter for simulated 
clusters at z = as functions of limiting mass. See the text for the 
definition of error bars. 
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Fig. 4. — Skewness and kurtosis of the Y—M scatter for simulated 
clusters with M500 ^ 2 X 10 13 Mq as functions of redshift. See the 
text for the definition of error bars. 



along the line of sight. Since less massive clusters are 
more likely to be surrounded by gas with mass compara- 
ble to their own, this effect becomes more important for 
lower-mass clusters. We will further address this point in 
§ 14.31 For other redshifts, the level of non-lognormality 
is also non- negligible, as shown in Figure |4j The median 
skewness and kurtosis are 1.43 and 4.21, respectively. 

4. SOURCES OF SCATTER 

We now investigate the physical origin of the intrin- 
sic scatter for the purpose of understanding the above 
trends and reducing it for better mass estimates. Pos- 
sible sources of scatter in our simulation include halo 
concentration, dynamical state, and cluster morphology. 
We examine each effect by correlating the scatter with 
each individual source for clusters at z = 0. We show 



that the scatter can be reduced by choosing appropriate 
measures of each effect. At the end of this section we 
compare the percentage contribution from each source 
using the idealized cluster samples described in § 12.41 

4.1. Concentration 

The concentration of a dark matter halo is usually 
defined as the ratio between the virial radius and the 
NFW scale radius, i.e., c = R V n/R s . The concentra- 
tion parameter characterizes the density inside the core 
region of a halo and reflects the mean density of the uni- 
verse when the halo collapsed. Thus halos forme d ear- 
lier in time ten d to be more concentrated (NFW Il997t 
iWechsler et all 2002). By testing the correlation of Y-M 
scatter with scatter in the concentration parameter, we 
are effectively probing the influence of cluster formation 
history. 

We choose to use the parameter -R200/-R500 instead of 
the original halo concentration parameter c because it 
has two advantages. The first is that it avoids introduc- 
ing the uncertainty of fitting an NFW profile, especially 
for less massive clusters, since the fitting is very sensitive 
to the grid resolution in the central region of a cluster. 
Moreover, our analyses involve not only relaxed clusters 
but also merging ones, for which R200/R500 is actually 
better-defined than c, since an NFW profile would yield 
a poor fit. 

Figure[5](lcft panel) shows a strong positive correlation 
between scatter in the Y-M relation and scatter in the 
(-R20o/-R50o) _ -^500 relation. The correlation coefficient is 
0.64, with a probability of zer o given by t he Spearman 
Rank-Order Correlation test (jPress et al.l 1992, §14.6; 
probability of one means no correlation) . To ensure that 
this result is not biased by the lower-mass clusters whose 
-R500 values are close to the resolution of the simulation, 
we raised the mass threshold to M 500 ^ 1O 14 M and 
found that the result is robust for these well-resolved 
systems (shown as the open squares in Figure [5]) . Note 
that we correlate with 51og(i?2oo/^50o) instead of the 
raw value of -R200/-R500 because the latter is a function 
of cluster mass. By doing so we exclude the effect of dif- 
ferent cluster masses, focusing on the variation in halo 
concentrations. -R200/-R500 is a monotonically decreas- 
ing function of t he halo concentration parameter (see 
lYang et all (|2009D for derivation). Therefore, for clusters 
with similar masses, more concentrated clusters tend to 
lie under the mean Y-M relation, while the "puffier" 
clusters tend to scatter high. 

Since halo concentration is related to cluster forma- 
tion time, we can test the above trend by checking the 
correlation of Y-M scatter with the formation times de- 
rived from the mass assembly histories of our simulated 
clusters. The formation time here is defined as the time 
when the cluster first exceeds half of its final mass. As 
expected, we find that clusters that formed earlier (thus 
with higher concentrations) tend to scatter low (see right 
panel of Figure [5]). The correlation is not as tight as the 
one with the halo concentrations. This is due to the fact 
that the correlation between the halo concentration and 
the cluster formation time itself has a very large scat- 
ter, and also that the variation in halo concentrations 
cannot be fully accounted for by the variation in cluster 
formation time ()Neto et al.ll2007t ). But the direction of 
the correlation with cluster formation time is consistent 
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Fig. 5. — Left: Y-M scatter versus (i?20o/^50o)~A^ scatter at z = 0, where -R20o/^?500 ls a monotonically decreasing function of halo 
concentration. Correlation coefficient is 0.775. Right: Y—M scatter versus formation time at z = 0. Correlation coefficient is -0.244. 
Clusters with M500 1O 14 M0 are plotted using open squares. 



with the correlation with halo concentration. 

To explain the correlation between the Y-M scatter 
and the concentration, recall the virial theorem for the 
simplest case of an isolated system: 2K + {7 = 0, where 
K and U are the total kinetic and gravitational binding 
energies of the system, respectively. In general one can 
write- 

k B T GM 

<x —5-, (12) 

/im p R 

where (i is the mean molecular mass of the gas, m p is the 
mass of a proton, and T , M, and R are the virial tem- 
perature, mass, and radius of the system, respectively. 
Together with the definitions of the SZ flux (Eq. [2]) and 
M = ^irR 3 p, one can derive 

Y cx M gas T cx f g M 5 / 3 . (13) 

Note that the above relations arc for virial quantities 
of a cluster as a whole, but mass-observable relations 
are often measured using a certain aperture size, Ra- 
The relation between the virial and overdensity quan- 
tities depends on individual cluster profiles, which are 
determined by the halo concentration and how the gas 
is distributed on top of the dark matter potential (e.g. 
equation of state of gas). Therefore, the normalization, 
and thus the scatter, of the Y&-M& relation is a function 
of halo concentration and gas properties. 

It is also important to note that, from the above deriva- 
tion, the direction of the correlation between the scatter 
and concentration is dependent on the gas properties. 
Our result shows that less concentrated clusters tend to 
scatter high, i.e., have higher p ressure than clusters o f 
similar masses. Observationallv. [Comerford et al.l (|2010f ) 
has also found a similar anti-correlation between the X- 
ray temperature-mass scatter and strong lensing concen- 
tration. This may be attributed to the fact that less 
concentrated clusters have larger scale radii, and hence 
when comparing with clusters of the same Ma or the 
same aperture size Ra, their ratios Ra/R s are smaller, 
which means the observable integrated within Ra would 
be greater. However, different simulations can have dif- 
ferent directions of correlation depending on the input 



gas physics. For example, iShaw et all (j2008| ) also found 
a correlation between the scatter and concentration, but 
in the opposite direction. The difference may be due to 
different gas physics included in their models. In prin- 
ciple, by assuming a particular gas model the constant 
of proportionali t y in Eq. 1131 can be computed exactly. 
lAscasibar et al.1 1)20061 ) has done this exercise assuming 
a polytropic equation of state for the gas. According to 
their calculation, the coefficient in the T-M relation (as 
in Eq. [TJ] inverse of the Ymt m their Eq. 17) decreases 
with concentration for a polytropic index of 7 P = 5/3. As 
7 P decreases, the dependence becomes weaker and then 
the direction is reversed. Since including extra bary- 
onic physics effectively works to decrease j p (e.g. for 
fixed mass and concentration, both changes yield a shal- 
lower tempera t ure p rofile, see also Figures 2 and 3 in 
lOstriker et all (120051 )). this may explain why the depen- 
dence on concentration can be different between models 
with different input gas physics. Note, however, that in 
reality the situation can be even more complicated be- 
cause a constant j p may not be valid for all gas models 
(IKav et al.ll2"00l . 

The strong correlation in Figure [5] suggests that the 
variation in halo concentrations contributes a signifi- 
cant amount of the intrinsic scatter in the Y-M rela- 
tion. Using this strong correlation it is possible to ad- 
just for the dependence of SZ flux on cluster concen- 
trations. We use <51og(i?20o/-R50o) for each cluster to 
calculate its expected 5\ogY from the best-fit relation, 
(51ogF) ex p = 7.167 x S\og(R 20 o/R50o)- We then sub- 
tract this (5 log F) eX p from the measured SZ flux to ob- 
tain a corrected flux, 

(S\ogY) COII = SlogY - (<51ogY)ex P . (14) 

The Y-M relations before and after correcting for con- 
centration are shown in Figure [6] We find that after 
removing the effect of halo concentration, the RMS scat- 
ter decreases from 12.07% to 7.34% (i. e by 3 8.9%). This 
method was proposed by lYang et all ([2009D to tighten 
the X-ray temperature-mass relation and has been suc- 
cessfully applied to o bserved strong lensing clusters 
(|Comerford et al.l 120101 ) . In addition to strong lensing, 
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Fig. 6. — Left: The Y-M relation plotted for clusters having different values of concentration at z = 0. The 1/3 percentiles with the 
highest, intermediate, and lowest values of concentration are plotted usin g black, blue, and red symbols, respectively. Right: The Y-M 
relation after correction for the dependence on concentration using Eq. 1141 



the NFW concentration can al so be measured via weak 
lensing, X-ray emission, etc. (IComerford fc Nataraianl 
120071: iMandeibaum et all l2008t iBuote et alj 120071 and 
references therein) , although one has to be careful about 
systematics of each method and when combining differ- 
ent measurements. In fact it can be generalized to any 
observable other than concentration. That is, if there ex- 
ists any variable X whose mass scatter S log X is known 
and correlates with S\ogY, then given the best- fit corre- 
lation (<5 log Y) cxp = ax 5 log X, the Y-M scatter can be 
reduced in a similar way using Eq. [14] to remove the ef- 
fect of X from the scatter. Therefore, this method can be 
a powerful way to reduce the observed mass-observable 
scatter and obtain better mass estimates. 

4.2. Dynamical State 

Another possible origin of the Y-M scatter is clus- 
ter dynamical state. Cluster mergers are among the 
most energetic events in the universe. Shock heating 
and departure from hydrostatic equilibrium during merg- 



ers can drive clusters away fro m the mean scaling rela- 
tions. Ideal merg er simulations (iRicker fc Sarazinll2001l : 
iPoole et al.l 120071 ) have shown that the effect of shock 
heating to boost the SZ and X-ray observables to values 
a few times higher than the pre-merger values. Stud- 
ies that combine the amount of boosting predicted by 
these simulations with extended Press-Sch echter merger 
trees ([Randall et al.ll2002tlWik et alJl2008| ) show that the 
boosting effect can bias estimates of cosmological param- 
eters such as <78 an d The other effect of mergers 
is departure from hydrostatic equilibrium. Before the 
gas within a merger is completely virialized, the pres- 
sure support from random gas motio ns can contribute 
~ 10 — 20% of its thermal pressure (|Rasia et al.l [2001 
lLau et al. 2009]). Thus the thermal pressure and hence 
the SZ flux of unrclaxed systems is expected to be smaller 
than that of relaxed systems of similar masses. These 
previous studies are primarily based on small cluster 
samples. Therefore, our aim is to investigate how merger 
events statistically influence the cluster scaling relations. 



If the scatter were dominated by the boosting effect of 
cluster mergers as described above, then one would ex- 
pect to find merging clusters to preferentially lie above 
the mean relation. However, if during mergers the de- 
parture from hydrostatic equilibrium due to non-thermal 
pressure support were dominant, mergers would tend to 
scatter low. In order to see which effect is more promi- 
nent, we correlate the Y-M scatter with the time since 
last merger (Figure [71 le ft panel). Substru cture measures 
such as centroid offset ([Mohr et aTlll995t) and power ra- 
tios (jBuote fc T sai 19951 11996| ) are often used to quantify 
departures from equilibrium in clusters. We compute the 
centroid offset and power ratios for the simu l ated c lusters 
using the same definition as in lYang et al.1 (pOOl . The 
right panel of Figure [7J shows the correlation with one of 
the power ratios, P2I Py- Based on th e Spearman Rank- 
Order Correlation test (jPress et al.l 1992, §14.6), both 
correlations are weak (with a small correlation coeffi- 
cient) but significant (with a high probability) , in the di- 
rection that more disturbed clusters tend to scatter low. 
This implies that during mergers the incomplete virial- 
ization may be the more important factor in driving the 
scatter than the shock heating effect. However, the fact 
that these two effects operate in opposite directions may 
be the reason why there is not a clear trend with merger 
activities. Moreover, a number of factors can dilute the 
shock boosting effect, such as the small chance of finding 
mergers in progress, capturing the shocks within an over- 
density radius at the right projection and at the right mo- 
ment during a merger's transient boost, and the fact that 
merging clusters tend to move along the scaling relations 
because their masses also increase at t he same time their 
obs ervables increase , as found also by IWik et alj pOOl 
and lKravtsov et al.l ([20061) (see lYang et al.1 (|2009f ) for an 
extensive discussion). 

Do mergers bias the Y-M relation due to incomplete 
virialization? To answer this question, we plot the nor- 
malized distributions of the Y-M scatter for relaxed and 
merging clusters in Figure [5] and use the Wilcoxon Rank- 
Sum (R-S) test and the F-variance (F-V) test to see 
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Fig. 7. — Left: Y—M scatter versus time since last merger for major mergers at z = 0. The correlation is statistically significant but 
weak (correlation coefficient of 0.299; probability of no correlation of 0.006). Right: Y-M scatter versus one of the power ratios, P2/ Pq- 
Correlation coefficient is -0.09f ; probability of no correlation is 0.023. 
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Fig. 8. — Normalized distribution of the Y-M scatter for relaxed 
(solid) and merging (dashed) clusters at 2 = 0. According to results 
from significance tests, mergers do not have a bias but do have a 
larger dispersion with respect to the relaxed clusters. 



whether these two populations have significantly different 
mean values or variances, respectively. A value smaller 
than 0.05 (for a significance level of 5%) returned by 
these tests is commonly adopted to indicate a significant 
difference between two populations. We find that their 
mean values do not differ significantly, but mergers have 
a wider distribution compared to relaxed clusters (with 
significance 0.014). Therefore, although mergers do not 
tend to bias the scaling relation, they do have a greater 
amount of scatter than relaxed clusters. If merging (re- 
laxed) clusters are chosen to be those within the quartile 
with the highest (lowest) substructure measures, similar 
trends with significant probabilities are found for 9 out 
of 21 substructure measures (P2/P0, P3/P0, and the cen- 
troid offset W measured from different viewing directions 
and varying aperture sizes). We find that including only 
the merging clusters would result in ~ 15 — 45% greater 
scatter than when only relaxed clusters ar e taken into 
accou nt, consistent with previous findings (|Shaw et al.l 
I2008T ) . Note however that separating mergers from re- 
laxed clusters docs not reduce the skewness or kurtosis 
of the scatter distribution, which suggests that th e no n- 
lognormality has causes other than mergers (see § I4.3j) . 



4.3. Morphology 

Despite the spherical symmetry that theo- 
retical models usually assume, both simulations 
dWarren et al.l fl992l: ICole fc Lacevi fT996t l.ling fc Sutol 
I2002t IBailin fc Steinmetzl 120051: IKasun fc Evrardl (20051) 
and observations (e.g. iBasilakos et al.ll2000[ ) have shown 
that clusters are triaxial (or elliptical when projected) 
rather than simple spheres, even for relaxed clusters. 
How the gas is distributed in the cluster potential well 
should vary depending on the axes ratios of the cluster. 
Moreover, viewing a triaxial cluster from different angles 
should also yield different observed quantities integrated 
along the line of sight. Both these factors can contribute 
to the Y-M scatter. 

In order to explore the impact of morphology, for each 
simulated cluster we find the orientation of the principal 
axes by diagonalizing the moment of inertia tensor, I a p = 

rriiYiirfr^ where the summation is over all the particles 
and cells in the cluster, uii is the mass of a particle or 
a gas cell, and rf is the x, y, or z component of the 
distance from the cluster center of mass. The lengths of 
the major, intermediate, and minor axes, denoted as a, 
6, and c, are found by finding the intercepts of the axes 
with the isodensity surface having overdensity A = 200. 
The angles 6 a are defined to be the angles between the 
major axis and the directions of projection (a = x, y, z; 
note that the x-direction is the projection used for all the 
analyses in this paper). 

Motivated by our results in § 13.31 that the non- 
lognormality may be due to clusters that happen to 
have elongated shapes aligned with the viewing direc- 
tion, we invented a measure, acos^, to trace cluster 
morphology along the line of sight. Figure O shows the 
Y-M scatter versus the scatter in the a cos 9 X -M re- 
lation. The positive correlation indicates that clusters 
that are more elongated along the line of sight prefer- 
entially have higher Y-M scatter, which is expected be- 
cause the SZ flux is roughly proportional to the column 
density of the gas (see Eq.[TJ. Given the best- fit relation, 
5 logy = 0.076 x 5 log (a cos 9 X ), we can again reduce the 
scatter by applying a correction as in Eq. 1141 By doing 
so we find that the scatter is reduced from 12.07% to 
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Fig. 9. — Correlation between the Y—M scatter and the acos9 x - 
M scatter at z = 0. Clusters that are more elongated along the 
line-of-sight have larger values of a cos 9 X . Correlation coefficient 
is 0.345. 
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Fig. 10. — Normalized distributions of the Y—M scatter for clus- 
ters with the morphology measure a cos X smaller than the me- 
dian (solid) and larger than the median (dashed). The scatter dis- 
tribution for clusters with elongated shape along the line-of-sight 
(dashed; 7 = 1.04, k = 2.63) is much more non-lognormal than 
that of the remaining population (solid; 7 = 0.29, k = —0.08). 

11.63% (i.e. by 3.6%). 

We further divide the cluster sample in half using the 
values of a cos 9 X and plot the normalized distributions of 
the Y—M scatter in Figure [TO] We find that the scatter 
distribution of the clusters with larger a cos 9 X is non- 
lognormal (7 = 1.04, k = 2.63), while the skewness and 
kurtosis of the remaining population are greatly reduced 
(7 = 0.29, k = —0.08). Therefore, the non-lognormality 
is indeed caused by the clusters with more elongated 
shapes along the line of sight. 

4.4. Projection Effects Due to Large-Scale Structure 

Since the SZ flux is obtained by integrating along the 
line of sight within a projected radius, it is subject to con- 
tamination by gas that lies along the same line of sight, 
which causes the large number of high-scatter objects in 
the Y- M relation found in simulations that incl ude light 
cones (jWhite et al.l [2001 lHallman et all [2007h . These 
outliers and the outliers due to morphology discussed in 
the previous section can both drive the non-lognormality 
of the scatter. Since our simulated observations only 




Fig. 11. — The percentage contribution of the scatter for each ide- 
alized sample with respect to the simulated sample. The bottom 
(top) sample includes the least (most) physical sources of varia- 
tions. See Table [T] for a summary of notations and assumptions 
used to construct each sample. 

include isolated clusters and thus do not take the pro- 
jection effect into account, the skewness and kurtosis es- 
timated from our simulation may be considered to be 
underestimates of the true values. In this case, it is 
even more important to adopt the higher moments in the 
parametrization of the scatter in order to get unbiased 
cosmological constraints. 

4.5. Insights from Idealized Samples 

Figure fTTI shows the percentage contribution of scatter 
for each idealized cluster sample with respect to the sim- 
ulated sample. The cluster sample at the bottom has the 
most constraints on and least freedom in the model pa- 
rameters, and the assumptions are loosened one at a time 
from bottom to top (see Table Q] for a summary of model 
descriptions). In general the values are independent of 
mass, that is, the processes shape the scaling relation in a 
self-similar way, as expected in the absence of additional 
baryonic physics. The only exception is the bottom curve 
for which only the masses of clusters are assigned. In 
principle this sample should have zero scatter if the res- 
olution of gas cells were infinite, but in reality the finite 
resolution introduces a nonzero scatter which becomes 
bigger as more low-mass clusters are included. Because 
the idealized samples are all constructed under the as- 
sumption of spherical symmetry, we derived Sample SS 
(second line from the top in the figure) using gas profiles 
directly extracted from the simulated clusters, such that 
it includes all sources of scatter except the morpholog- 
ical effect. In other words, the difference between the 
simulated clusters and the spherically-smoothed clusters 
is solely due to the variation in cluster morphology, which 
is - 10%. 

From the differences between the subsequent samples 
we are able to isolate the contribution of each effect 
to the total scatter: the variation in halo concentra- 
tion contributes ^ 10 — 20 % (difference between A and 
B), the departure from hydrostatic equilibrium results 
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in ~ 10 — 15% (between B and C), merger boosts add 
another ~ 30 — 60 % (between C and D), the varia- 
tion in gas fractions introduces ~ — 10% (between D 
and E), and the rest (between D and SS) due to other 
unaccounted-for effects is ~ — 30 %. Note that the con- 
tribution from variation in concentration quoted here is 
estimated using sample B, which assumes spherical sym- 
metry and hydrostatic equilibrium. However, in reality, 
both changin g cluster morpho logy and including random 
gas motions (jLau e t al. 20091) would further modify the 
distribution of gas and hence alter the measured value of 
concentration. In other words, the scatter is driven not 
only by the variation of concentration in sample B, but 
also by that due to cluster morphology and departure 
from HSE. Therefore, the total effect of concentration, 
as suggested by the correction for concentration in § 14.11 
(i.e. ~ 40%), would be more appropriately accounted for 
by also considering the contributions from morphology 
and random gas motions. 

5. COMBINING X-RAY AND SZ SCALING RELATIONS 

In the previous sections we have discussed various 
sources of intrinsic scatter in the relation between the SZ 
flux and the true mass. However, observationally cluster 
masses still need to be measured in some way, such as via 
X-ray hydrostatic assumptions or optical richness. Cross- 
calibration across measurements at different wavelengths 
is important because it provides a consistency check that 
can minimize the possible s ystematic effects o f each in- 
dividual measurement (e.g. IPlagge et al.ll2010h . such as 
the projection effects to which SZ and optical observa- 
tions are subjected to. Therefore, high-precision clus- 
ter cosmology requires that we combi ne SZ cluster sur- 
veys with X-ray or op tical follow-ups (jffigh et fdi&M 
iMenanteau et alJl2010h . 

However, one needs to be cautious when combining 
multiple mass proxies because their errors may be cor- 
related. For example, because both the SZ and opti- 
cal signals arc subject to projection effects, clusters can 
have consistent mass estimates th at are both actually b i- 
ascd with respect to the true mass (jCohn fc Whitd[2009h . 
Since it is impossible for observations to disentangle such 
correlations, one has to rely on numerical simulations 
to determine whether these effects are serious for any 
given pair of mass proxies. Here we would like to explore 
whether this correlated error exists betw een the SZ flux 
and the low-scatter X-ray mass proxy (|Kravtsov et al.l 
|2006[ ). Yx, which is commonly used as a mass proxy in 
X-ray observations. Note that although individual X-ray 
properties such as M gas and Tx would be affected by core 
properties, Stanek et al. (2007) found that the Yx pa- 
rameter, which combines the effects of M gas and Tx, is 
remarkably insensitive to baryonic physics. That is, for 
their runs with and without the preheating prescription, 
both the amount and shape of the Yx~M scatter are al- 
most identical. Therefore, our results below should be 
robust to additional baryonic physics. 

Figure [T2l shows the mass predicted by the Y$z~M re- 
lation versus that predicted by the Yx~M relation (M 
is the true mass). Clusters in the upper panel have less 
than 2cr deviations from both the mean Ysz~M relation 
and Yx~M relation. The lower panel shows the clusters 
whose mass scatter is bigger than 2a for cither relation. 
From the upper panel we can see that clusters that have 
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Fig. 12. — Ysz predicted mass versus Yx predicted mass. Clus- 
ters in the upper panel have less than 2<r deviations from both the 
mean Yg^—M relation and Yx~ M relation, while clusters whose 
mass scatter is bigger than 2<r for either relation are plotted in 
the lower panel. Dashed lines show lcr deviations from the mean 
My sz —My x relation. The fact that clusters that are outliers in 
both relations (those with overlaying circle and triangle) do not 
have consistent mass estimates within lcr indicates that the errors 
in My sz and My x are not correlated. 



consistent My sz and My x are mostly faithful tracers of 
their true masses. But how about the outliers in both the 
true Ysz-M and Yx~M relations? If they give consistent 
mass estimates, then there would be a similar problem of 
correlated error as described above. Fortunately, we find 
that the outliers in both relations (those plotted with 
both open and filled symbols) would not yield consis- 
tent mass estimates. This is because while the Ysz~M 
outliers are due to cluster morphology, as discussed in 
§ I4.3[ we nn d that the Yx~M outliers are primarily dy- 
namically unrelaxed clusters. Since the errors come from 
different physical sources, they are not correlated. 

This implies the possibility of cutting off the outliers 
by selecting only clusters whose My sz and My x agree 
within lcr. Moreover, applying the same cut will also 
remove almost all the other Y$z~M outliers. That is, 
among the 21 Ysz~M outliers (15 are 2a and 6 are 3cr), 
19 of them (13 are 2a and 6 are 3a) can be ruled out 
using this method. After applying the cut, we find that 
the RMS scatter in Y SZ ~M is reduced from 12.07% to 
8.77% (i.e. by 27.3%), and also the non-lognormality of 
the Ysz-M scatter is greatly reduced (skewness reduced 
from 0.82 to 0.30; kurtosis from 2.17 to -0.07). There- 
fore, combining mass estimates from Yx measurements 
may be an effective way of both reducing the scatter and 
removing Ysz~M outliers. Note that since the projec- 
tion effect is not included in our simulation, we expect 
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there would be more Ysz~M outliers in reality, while the 
Yx~M relation is relatively insensitive to the projections. 
In deed, the contaminatio n by projection errors estimated 
by lHallman et all (|2007j ) using a light cone simulation is 
~ 25% (for a projected radius of -R500), larger than ours 
(21 out of 619 clusters). However, because of the fact 
that the errors in the Ysz~M and Yx~ M relations are 
not correlated, clusters that are subject to projection er- 
rors can also be removed using the same method. 

6. DISCUSSION AND CONCLUSIONS 

Galaxy clusters are invaluable cosmological probes. 
Accurate measurement of cluster masses is crucial and of- 
ten relies on the mass-observable relations. However, to 
constrain the cosmological parameters at the few percent 
level, the systematics and scatter in these relations must 
be thoroughly understood. In this work we investigated 
the sources of intrinsic scatter in the SZ flux-mass (Y-M) 
relation using a hydrodynamics plus iV-body simulation 
of galaxy clusters within a cosmological volume. Explor- 
ing the origin of the intrinsic scatter not only provides 
physical insights into the formation of galaxy clusters, 
but also has two main advantages for using clusters in 
cosmology. The first is that it allows us to avoid pos- 
sible systematic biases in the derived cosmological con- 
straints. Do mergers bias the scaling relation? Is the in- 
trinsic scatter lognormal? What are the gains and issues 
of combining SZ and X-ray scaling relations? Secondly, 
if we understand the sources of scatter, it is possible to 
reduce the scatter by removing the contribution from a 
certain source (see § 14.11 for details) , and thus tighten 
the scaling relation to obtain better estimates of cluster 
masses. 

To address these questions, we derived the scatter 
around the best-fit Y-M relation from the simulated 
clusters. We first assessed the lognormality of scatter 
by computing the skewness (7) and kurtosis (k) of the 
scatter distribution. Then we investigated the possible 
sources of scatter, including halo concentrations, dynam- 
ical states and cluster morphology, by correlating the 
scatter with quantitative measures of each source. We 
also constructed a set of idealized cluster samples with 
varied assumptions about the sources of scatter to de- 
compose the percentage contribution from each effect. 
Finally we compared cluster masses derived from the SZ 
flux and from the low-scatter X-ray mass proxy, Yx , and 
examined whether such consistency checks can help rule 
out outliers in the true Ygz~M and Yx~M relations, or 
whether issues like correlated errors would affect the ac- 
curacy when combining SZ and X-ray scaling relations. 
Our main results are summarized below. 

1. The RMS scatter in the Y-M relation is ~ 5 - 15% 
and decreases with cluster masses and redshifts. We find 
that the scatter in our simulation can be expressed in the 
functional form, a(M,z) = AlogM + Blog(l + z) + C 
(Eq. [7]), where the redshift evolution is equivalent to re- 
scaling with respect to the characteristic mass scale in 
the self-similar model. 

2. The distribution of the Y-M scatter is non- 
lognormal with positive skewness and kurtosis across a 
wide range of different limiting masses and redshifts, be- 
cause of the limited number of clusters at the higher- mass 
end and the tail in the scatter distribution due to mor- 
phology at the lower-mass end. 



3. There is a strong correlation between the Y-M scat- 
ter and the concentration, which can be used to reduce 
the Y-M scatter from 12.07% to 7.34% (i.e. by 38.9%). 

4. The correlation between the scatter and cluster dy- 
namical state is weak. Though merger boosts and depar- 
ture from hydrostatic equilibrium can partly drive the 
dispersion, the net effect is that mergers do not cause a 
significant bias in the scaling relation. 

5. There is a moderate trend that clusters that are 
more elongated along the line of sight tend to scatter 
high. More importantly, they are the main outliers that 
cause the non-lognormality of scatter. 

6. By decomposing the scatter using the idealized clus- 
ter samples, we find the percentage contribution from 
each source of scatter: ~ 10% due to variations in mor- 
phology, ~ 10 — 20% due to variations in concentration 
(under the assumption of spherical symmetry and hydro- 
static equilibrium), ~ 10—15% due to departure from hy- 
drostatic equilibrium, ~ 30 — 60% due to merger boosts, 
~ — 10% from variations in gas fractions. The remain- 
der (due to unaccounted-for sources) is ~ — 30%. 

7. We find that the RMS scatter in Y sz -M is reduced 
from 12.07% to 8.77% (i.e. by 27.3%) when X-ray mea- 
surements are combined with SZ. 

8. The errors in mass determined using Ygz arL d Yx 
come from different causes. Therefore, excluding clus- 
ters with inconsistent estimates can effectively remove 
the outliers in both Ysz~M and Yx~M relations, es- 
pecially Ysz-M outliers that are subject to projection 
errors. 

In our current simulation, radiative cooling and heat- 
ing mechanisms are not included, since we would like 
to disentangle the scatter driven by the gravitational ef- 
fects from other baryonic physics that are not fully un- 
derstood. Moreover, it has been shown that the inte- 
grated SZ flux, and more specifically the scatter, slope, 
and redshift evolution of the Y-M relatio n, are generally 
insen s itive to details of c l uster gas phy sics ()da Silva et al.l 
l200i IMotl et all [20051: INaeail 120061) . Since the non- 
lognormality is mainly caused by the effects of projec- 
tion and cluster morphology, we could assess the po- 
tential impact of baryonic physics on these two sources. 
iShaw et al.l (|2008t ) showed that the influence of different 
gas physics on the pr operties of large -scale projections 
is negligible. Recently lLau et ail (|2010f) has reported the 
difference in cluster shapes between simulations with and 
without cooling and star formation. They found that in 
the cooling plus star formation simulation clusters arc 
more spherical outside the core (r > O.I-R500) but more 
triaxial inside the core. It is difficult to estimate directly 
from their results how much this difference in morphol- 
ogy would affect the non-lognor mality of the Y-M scat- 
ter. However, as pointed out by ILau et al.l (|2010D . their 
simulation may suffer from the overcooling problem and 
hence thei r resul t s can be considered as an upper limit. 
Moreover . iNagail (|2006t ) used clusters from the same sim- 
ulations and showed that when the SZ flux is integrated 
to -R500 , the Y-M scatter is insensitive to the gas physics 
included. Therefore, we expect the effect of gas physics 
on the non-lognormality estimated in this paper, if any, is 
very small. We will present a more detailed comparison 
in a separate paper. 

Our results have several important implications for 
cluster cosmology. First of all, the strong correlation 
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with halo concentrations can be used for observed clus- 
ters to reduce the scatter in the scaling relations for 
better mass estimates. Potentially this method can be 
applied to any observable for which such a correlation 
exists, such as gas fractions (which is expected to play 
a more imp ortant role when in cluding other baryonic 
physics, see iStanek et al.l (|2010D ). Secondly, the weak 
influence of mergers is good news for cluster cosmology, 
because it implies that when deriving observed scaling 
relations, it is unnecessary to worry much about the 
selection bias due to the impact of mergers. Finally, 
the non-lognormality of the Y—M scatter has an impact 
on cosmological con straint studies. As demonstrated by 
iShaw et ajX(|2010al) . both positive skewness and kurto- 
sis cause up-scattering of clusters and thus would in- 
crease cluster counts above a given limiting mass, which 
is equivalent to an increase in the amount of scatter. 
For SZ surveys like the South Pole Telescope (SPT) sur- 
vey, the skewness and kurtosis of the intrinsic scatter 
have to be less than 0.5 to ensure that uncertainty in the 
amount of scatter does not degrade the constraint on the 
dark energy equation of state w. However, we find that 
the intrinsic skewness and kurtosis can be much greater 
than 0.5 across a wide range of limiting masses and red- 
shifts. These values are very likely to be lower limits be- 
cause the projection effect of large-scale structure is ab- 
sent in our analysis. Therefore, our results suggest that 
the assumption of lognormal scatter is inappropriate for 



scaling relations like the Y—M relation whose scatter is 
easily skewed by cluster morphology, projection effects, 
etc. Instead, in self-calibration studies of cosmological 
constraints that require assumptions about the form of 
scatter, it is necessary to include the higher-order mo- 
ments in the parametrization. During the next decade, 
more and more data from multi- wavelength cluster obser- 
vations will be available. We expect that more detailed 
studies of the intrinsic scatter in the scaling relations will 
continue to yield essential information both for cluster 
physics and cluster cosmology. 
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